# Goal: Make Table D1
# Dependencies: "datamodel.Rdata", "reffs_samples.Rdata"

setwd("~/Replication Package")

library(coda)
library(gplots)
library(xtable)

# Load survey data
load("datamodel.Rdata")

datamodel2 <- data.frame("education" = datamodel$education, "age" = datamodel$age, "female" = datamodel$female, "nonwhite" = datamodel$nonwhite, "income.proxy" = datamodel$income.proxy, "ideology" = datamodel$ideology, "natecon" = datamodel$natecon, "persfin" = datamodel$persfin, "victim" = datamodel$victim, "corrupt" = datamodel$corrupt, "capital" = datamodel$capital, "bigcity" = datamodel$bigcity)

# Standardize covariates
datamodel2 <- as.data.frame(scale(datamodel2))

# Load 2nd stage chains
load("reffs_samples.Rdata")

reffs.samples <- as.mcmc.list(reffs.samples)

coef.samples <- reffs.samples[ , substr(names(reffs.samples[[1]][1, ]), 1, 12) != "Sigma.beta[1" & substr(names(reffs.samples[[1]][1, ]), 1, 6) != "beta[1" & substr(names(reffs.samples[[1]][1, ]), 1, 9) != "mu.beta[1"]

coef.samples.all30 <- NULL
for (j in 1:30) {
  coef.samples.all30 <- rbind(coef.samples.all30, coef.samples[[j]])
}

N.iters <- dim(coef.samples.all30)[1]

n.sims <- 13

X.median <- c(1, apply(datamodel2, 2, median))
names(X.median)[1] <- "constant"
X.median["female"] <- quantile(datamodel2$female)[2]
X.median["natecon"] <- quantile(datamodel2$natecon)[5]
X.median["persfin"] <- quantile(datamodel2$persfin)[5] 
X.median["corrupt"] <- quantile(datamodel2$corrupt, p = c(0.05)) 

X.sims <- matrix(rep(X.median), nrow = n.sims, ncol = length(X.median), byrow = T)
colnames(X.sims) <- names(X.median)

X.sims[2, "education"] <- quantile(datamodel2$education)[4]
X.sims[3, "age"] <- quantile(datamodel2$age)[4]
X.sims[4, "female"] <- quantile(datamodel2$female)[4]
X.sims[5, "nonwhite"] <- quantile(datamodel2$nonwhite)[4]
X.sims[6, "income.proxy"] <- quantile(datamodel2$income.proxy)[4]
X.sims[7, "ideology"] <- quantile(datamodel2$ideology)[4]
X.sims[8, "natecon"] <- quantile(datamodel2$natecon)[2]
X.sims[9, "persfin"] <- quantile(datamodel2$persfin)[2]
X.sims[10, "victim"] <- quantile(datamodel2$victim)[5]
X.sims[11, "corrupt"] <- quantile(datamodel2$corrupt)[4]
X.sims[12, "capital"] <- quantile(datamodel2$capital)[5]
X.sims[13, "bigcity"] <- quantile(datamodel2$bigcity)[5]

betas.list <- list()
eXB.list <- list()
P.list <- list()

for (j in 2:4) {
  
  betas.list[[j]] <- array(coef.samples.all30[, substr(colnames(coef.samples.all30), 1, 6) == paste("beta[", j, sep="")], c(N.iters, 3, 13))
  
  eXB.list[[j]] <- array(NA, c(N.iters, 3, 13))
  
  for (i in 1:3) {
    eXB.list[[j]][ , i, ] <- exp(betas.list[[j]][ , i, ] %*% t(X.sims))
  }
  
}

for (j in 1:4) {
  P.list[[j]] <- array(NA, c(N.iters, 3, 13))
}

for (i in 1:3) {
  P.list[[2]][ , i, ] <- eXB.list[[2]][ , i, ] / (matrix(1, nrow = N.iters, ncol = 13) + eXB.list[[2]][ , i, ] + eXB.list[[3]][ , i, ] + eXB.list[[4]][ , i, ])
  P.list[[3]][ , i, ] <- eXB.list[[3]][ , i, ] / (matrix(1, nrow = N.iters, ncol = 13) + eXB.list[[2]][ , i, ] + eXB.list[[3]][ , i, ] + eXB.list[[4]][ , i, ])
  P.list[[4]][ , i, ] <- eXB.list[[4]][ , i, ] / (matrix(1, nrow = N.iters, ncol = 13) + eXB.list[[2]][ , i, ] + eXB.list[[3]][ , i, ] + eXB.list[[4]][ , i, ])
  P.list[[1]][ , i, ] <- matrix(1, nrow = N.iters, ncol = 13) -  P.list[[2]][ , i, ] -  P.list[[3]][ , i, ] -  P.list[[4]][ , i, ]
}

Base.mean.list <- list()
Base.ql.list <- list()
Base.qu.list <- list()
dP.list <- list()
dP.mean.list <- list()
dP.ql.list <- list()
dP.qu.list <- list()

for (j in 1:4) {
  
  Base.mean.list[[j]] <- apply(P.list[[j]][ , , 1], 2, mean)
  
  Base.ql.list[[j]] <- apply(P.list[[j]][ , , 1], 2, quantile, p = 0.025)
  
  Base.qu.list[[j]] <- apply(P.list[[j]][ , , 1], 2, quantile, p = 0.975)
  
  dP.list[[j]] <- array(NA, c(N.iters, 3, 12))
  
  for (k in 1:12) {
    dP.list[[j]][ , , k] <- P.list[[j]][ , , k+1] - P.list[[j]][ , , 1]
  }
  
  dP.mean.list[[j]] <- apply(dP.list[[j]] ,c(2, 3), mean)
  dP.ql.list[[j]] <- apply(dP.list[[j]] ,c(2, 3), quantile, p = 0.025)
  dP.qu.list[[j]] <- apply(dP.list[[j]] ,c(2, 3), quantile, p = 0.975)
  
  colnames(dP.mean.list[[j]]) <- names(X.median)[-1]
  colnames(dP.ql.list[[j]]) <- names(X.median)[-1]
  colnames(dP.qu.list[[j]]) <- names(X.median)[-1]
  
}

TableD1 <- list()

for (y in 1:3) {

  TableD1[[y]] <- cbind(
  dP.mean.list[[1]][y, ], dP.ql.list[[1]][y, ], dP.qu.list[[1]][y, ],
  dP.mean.list[[2]][y, ], dP.ql.list[[2]][y, ], dP.qu.list[[2]][y, ],
  dP.mean.list[[3]][y, ], dP.ql.list[[3]][y, ], dP.qu.list[[3]][y, ],
  dP.mean.list[[4]][y, ], dP.ql.list[[4]][y, ], dP.qu.list[[4]][y, ])
  
  TableD1[[y]] <- rbind(c(Base.mean.list[[1]][y], Base.ql.list[[1]][y], Base.qu.list[[1]][y], 
  Base.mean.list[[2]][y], Base.ql.list[[2]][y], Base.qu.list[[2]][y], 
  Base.mean.list[[3]][y], Base.ql.list[[3]][y], Base.qu.list[[3]][y], 
  Base.mean.list[[4]][y], Base.ql.list[[4]][y], Base.qu.list[[4]][y]), TableD1[[y]]) * 100

  colnames(TableD1[[y]]) <- c("P(Outsider)", "2.5%", "97.5%", "P(Agitator)", "2.5%", "97.5%", "P(Conventional)", "2.5%", "97.5%", "P(Activist)", "2.5%", "97.5%")
}

names(TableD1) <- c("2010", "2012", "2014")

xtable(TableD1[[1]], digits = 1) # 2010
xtable(TableD1[[2]], digits = 1) # 2012
xtable(TableD1[[3]], digits = 1) # 2014
